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1 Attorney Docket No. 83597 



2 

3 A METHOD FOR TRACKING TARGETS WITH HYPER- SPECTRAL DATA 

4 

5 STATEMENT OF GOVERNMENT INTEREST 

6 The invention described herein may be manufactured and used 



7 by or for the Government of the United States of America for 

8 governmental purposes without the payment of any royalties 

9 thereon or therefore . 
10 



11 CROSS REFERENCE TO OTHER PATENT APPLICATIONS 

12 Not applicable. 
13 

14 BACKGROUND OF THE INVENTION 

15 (1) Field of the Invention 

16 The present invention relates to remote sensing and remote 

17 imaging, and more specifically to a method for processing hyper- 

18 spectral remote sensor data for the purpose of displaying the 

19 spatial tracks of energy sources in multi-spectral images 
2 0 corresponding to the sensor data. 

21 (2) Description of the Prior Art 

22 Remote sensing of the energy signals of a moving vehicle or 
23^ energy source for the purposes of tracking the vehicle or energy 
24 source has often been accomplished by measuring the intensity of 
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1 the energy signals with sensors specifically designed to detect 

2 energy intensity. In remote sensing applications the sensor is 

3 often a planar array of sensing cells, each cell responding to 

4 the energy incident on its corresponding section of the array 

5 surface., In other applications, such as acoustic sensing, the 

6 received energy on sensor elements must be interpreted through a 

7 beam forming function to yield energy intensity in a set of 

8 spatially directed cells (more commonly called beams) . Such a 

9 method is designed to track energy peaks as they move over time 

10 on the given set of sensor cells. The total broadband energy is 

11 plotted and visually displayed. Targets appear as peaks of 

12 energy in the display, and are tracked. One method of target 

13 tracking based on sensor level data is the Histogram 

14 Probabilistic Mult i -Hypothesis ' Tracking (H-PMHT) algorithm. It 

15 is an application of the Expectation-Maximization (EM) method of 

16 target tracking. It uses a synthetic (multi-dimensional) 

17 histogram interpretation of the received power levels in all of 

18 the sensor cells. The data for the H-PMHT algorithm usually 

19 consists of broadband intensities on a set of spatial sensor 

2 0 cells. The H-PMHT algorithm has its limitations. For example, 

21 in situations where more than one target is being tracked and 

22 the targets cross paths, the intensity of the energy signals of 

23 the targets merge, making it impossible to distinguish the 

24 energy between the two targets. In such a situation, the 
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1 targets must be reacquired by the sensors after they have 

2 crossed resulting in a gap and delay in tracking information. 

3 U.S. Patent Application 10/214551 to Struzinski teaches a 

4 method and system for predicting and detecting the crossing of 

5 two target tracks in a bearing versus time coordinate frame. 

6 The method/system uses a series of periodic bearing measurements 

7 of the two target tracks to determine a bearing rate and a 

8 projected intercept with a bearing axis of the bearing versus 

9 time coordinate frame. A crossing time t c for the two target 

10 tracks is determined using the tracks 1 bearing rates and 

11 projected intercepts. A prediction that the two target tracks 

12 will cross results if a first inequality is satisfied while a 

13 detection that the two target tracks have crossed results if a 

14 second inequality is satisfied. This method does not, however, 

15 address the problem of distinguishing between and identifying 

16 the targets before, during and after they have crossed. 

17 There is currently no reliable method by which targets can 

18 be consistently tracked and distinguished as they cross paths. 

19 What is needed is a method for tracking targets that does not 

20 rely solely on broadband energy signal intensity, but also 

21 utilizes the spectral aspects of the energy signal, combining 

22 both intensity and spectral data so that crossing targets can be 

23 tracked provided they have some degree of spectral distinction. 
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1 SUMMARY OF THE INVENTION 

2 It is a general purpose and object of the present invention 

3 to provide a method for tracking both the spatial sensor data 

4 and hyper- spectral sensor data associated with a target. 

5 It is a further purpose to estimate a frequency spectrum 

6 for a target contribution that consists of only the target's 

7 energy contribution as opposed to the target energy and the 

8 noise energy together. 

9 These objects are accomplished with the present invention 

10 by taking the histogram model used in H-PMHT and extending it to 

11 treat the problem of tracking using hyper-spectral data. In the 

12 present invention each measurement scan is now a multi- 

13 dimensional array wherein each spatial cell has an associated 

14 vector of amplitudes in several (possibly disjoint) spectral 

15 cells. The intensity data in the multi-dimensional array is 

16 interpreted as the spatial -spectral histogram of a synthetic 

17 shot process. A statistical model of the random variation of 

18 individual cell intensities from scan to scan is required. The 

19 procedure adopted in H-PMHT is to quantize the data vector into 

20 a ''pseudo-histogram," and then use a multinomial distribution to 

21 model the cell counts where the PMHT target mixture model 

22 parameterizes the multinomial distribution. The target mixture 

23 model determines the cell probabilities that correspond to 

24 expected cell counts. 
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1 The present invention modifies H-PMHT by using a non- 

2 parametric spectral characterization of the energy intensity of 

3 the target that is assumed known. The use of such a spectral 

4 template enhances low signal to noise ratio (SNR) tracking and 

5 allows discrimination of spectrally distinct sources as they 

6 cross in the spatial domain. The track solutions from previous 

7 batches are used to estimate the (non-parametric) spectral 

8 characterization that is used to initiate the generation of the 

9 updated solutions as new data is received and processed. 

10 The present invention provides a mechanism for separating 

11 the observed hyperspectral energy into the hyperspectral energy 

12 for each source (including noise) using the known spectral 

13 characteristics for each source. Completely general spectral 

14 density functions are handled via the use of non-parametric 

15 methods. In the alternative, the source spectrum is estimated 

16 in a non-parametric fashion based on an initial track, allowing 

17 the algorithm to adapt to the source spectrum in situ. 
18 

19 BRIEF DESCRIPTION OF THE DRAWINGS 

2 0 A more complete understanding of the invention and many of 

21 the attendant advantages thereto will be readily appreciated as 

22 the same becomes better understood by reference to the following 

23 detailed description when considered in conjunction with the 
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1 accompanying drawings depicting an underwater application as the 

2 preferred embodiment wherein: 

3 FIG. 1 shows an underwater vehicle towing sensors; 

4 FIG. 2 shows sensors detecting energy from different 

5 sources; 

6 FIG. 3 shows the data cube created after raw sensor data is 

7 processed; and 

8 FIG. 4 shows a flow chart of the method. 
9 

10 DESCRIPTION OF THE PREFERRED EMBODIMENT 

11 Referring now to FIG. 1 there is shown an underwater 

12 vehicle 10 towing an array of sensors 15 arranged on a cable 20. 

13 The sensors 15 are of a type known by those skilled in the art 

14 of signal processing such as hydrophones. The sensors 15 are 

15 capable of detecting energy signals and their intensities from 

16 different directions as illustrated in FIG. 2, which shows two 

17 vessels labeled ki and k 2 and the energy signals 17 emanating 

18 from the vessels. The sensor data from each sensor 15 is 

19 transmitted along cable 20 to data processors (not shown) within 

20 underwater vehicle 10. The data processors take the raw sensor 

21 data and create a data cube 25 as illustrated in FIG. 3. Each 

22 such data cube 25 is a collection of smaller cubes referred to 

23 as sensor cells 30 which correspond to the processed sensor data 

24 generated by the sensors 15. Each sensor cell 30 contains 
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1 spatial measurements along the x-axis 31, spectral measurements 

2 along the y-axis 32 and time measurements along the t-axis 33. 

3 The side of each sensor cell 30 contained in the (x,t) plane 

4 corresponding to spatial measurement is referred to as a spatial 

5 cell 36. The side of each sensor cell 30 contained in the (y,t) 

6 plane corresponding to spectral measurement is referred to as a 

7 spectral cell 38. The processing and arrangement of the raw 

8 sensor data into a data cube 2 5 composed of multiple sensor 

9 cells 3 0 that are further composed of spatial cells 36 and 

10 spectral cells 38 is known by those skilled in the art of signal 

11 processing and is achieved by what is often termed beamforming 

12 followed by spectral analysis of the beam intensity data using 

13 standard discrete Fourier transform (DFT) techniques known in 

14 the art. A single layer of the data cube 25 is referred to as a 

15 scan of the sensor space 35 as illustrated in FIG. 3. 

16 In the preferred embodiment of the present invention, let 

17 C={C lf . . . , C s ), S > 1, denote the collection of all possible sensor 

18 cells 3 0.. It is assumed that C± n Cj = <t> for all i and j and 

19 that C x u.-.uCj = ^ din * c ) f W here dim(C) denotes the dimension of the 

20 sensor space 35. Furthermore, the sensor cells 30 C={Ci, . . 

21 .,C S } are the Cartesian products of U disjoint spatial cells 36 

22 {£>i. . and V disjoint spectral cells 38 {8i, . . ,,8 V }. This 

23 particular choice of spatial 36 and spectral 38 cells is 
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1 application dependent, but they are intrinsically fixed. The 

2 total number of sensor cells 3 0 in a scan of the sensor space 3 5 

3 is S=UV, and every cell Cj can be written in the form 

4 Cj = D ± xSj 

5 for some (unique) choice of the cells £>i and 6j .. Let D=£)iU. . 

6 .U £>u and 8=6i u. . . U8 V . The sensor cells 3 0 from which 

7 measurements are available may vary from scan to scan. The 

8 sensor cells 30 displayed at time t, for the scan of the sensor 

9 space 35, are the Cartesian product of the spatial cells 36 

10 {Di(t),-. . . ,D u(t ) (t) } and the spectral cells 38 {Ei(t),. . 

11 . ,E V (t) (t) } , so that the (i,j) th sensor cell 30 is Cij(t)=Di(t) x 

12 Ej(t), where 1 < U(t) < U and 1 < V(t) < V. The remaining sensor 

13 cells 30 in the scan of sensor space 35 are said to be 

14 truncated, and no measurements are collected for these cells at 

15 time t. 

16 Let the scan of sensor space 3 5 at time t be denoted by 

17 z t = \ Z t,l,l' • • • / z t,u{t) ,v(t)\ \ 

18 where z t ij > 0 is the output of the sensor space 3 5 at time t in 

19 cell CijCt), i = l,2,. . .U(t), j=l,2,\ . .V(t) . Let h 2 > 0 be a 

20 specified quantization level, and let n t ±j denote the quantized 

21 value corresponding to the intensity z ti j in cell Cij(t), where 



22 , ' n tij 



z tij 

h 2 



(i) 



1 and LxJ denotes the greatest integer less than or equal to x. 

2 The use of the quantized values {n tij } instead of the 

3 measurements {z ti j} is an intermediate step in the development. 

4 After deriving the auxiliary function of the H-PMHT algorithm 

5 using the synthetic counts {n t ij} , the measurements {z ti j} are 

6 recovered in the limit as h 2 — > 0 . 

7 The "rectangular" spatial -spectral sensor cell structure 

8 enables simplifications to the basic equations of H-PMHT. These 

9 equations are restated here with the updated notation 

10 corresponding to this new cell structure. The cell probability, 

11 Pij, for the (i,j) th cell 30 takes the form 

12 Pij&t) = I f(u f v\K t )du,dv f (2) 

13 where the sample Probability Density Function (PDF) f{u f v\x t ) is 

14 defined over all (u,v) e R dW? x R dim£ =R dimC , by the mixture 

15 density 

M 

16 f(u,v\x t )= J] Xtk G k&, v\x tk ) (3) 

k = o 

17 and where n tk is the component mixing proportion, 

18 X t = { x tO'*"' x tM} are the component spatial state parameters and 

19 G/c(u,v|x t k) is the component PDF corresponding to target 7c if k > 

20 1 and to noise if k = 0. The expected sensor space measurement 

21 Ztij takes the form 
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ztij 



ztij 



IN 



1 < i < U(t) , 

1 < j < V(t) , 

U(t) + 1 < i < U, 

V(t) + 1 < j <> V, 



(4) 



11 
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5 where |z t | = E^i z tiJ ' 



U(t) V{t) 

i=ij=i 



(5) 



7. and X t ' is the last estimate of X t . Thus, from (4) it may be 

8 seen that expected measurements exist for all cells, even those 

9 truncated in the observation. After taking the quantization 

10 limit, h 2 -> 0, the H-PMHT auxiliary functions become 



M 

Qt* = E 



t7 V 



12 and 



T 17 V 



^t* 109 ^t)c 



(6) 



(7) 



14 The density P~ i_ ( x tJc| x t-i,A) f° r t=l,2,. . .T describes the 



15 Markov process for the state of target k. 
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Let the spectral PDF of target k be denoted by iS k (v), so 



8 

9 
10 
11 

12 
13 
14 



that 



3 ^S k (y) dv = 1. (8) 

4 The spectral PDF is equal to the traditional power spectrum 

5 normalized so that its integral over e is one. Because the target 

6 spatial and spectral characteristics are independent by 

7 assumption, each component G k {u,v\yL tk ) of the sample PDF factors: 

G k (u, v\x tk ) = 0 k (u\x tk )S k (v) , (9) 



where g k (u\y: tk ) is the spatial PDF of component k. Independence 
enables integrals over Cij ( t) to be rewritten as products of 
integrals, so that 

l. jit) G k {u '^ x tk )du dv = L (t) ^ (v)dv t.( t) ^H x ^ )da (io) 



and, using the mixture (3) and the definition (2) , 

M 

p ij {x t) = z^i. (t) ^ (v)dv t. (t) ^H x ^ )da 

k = 0 J 1 



(id 



15 Substituting (10) into (6) gives 



M 



16 



17 where 



18 



k = 0 



2 V t ki\ i(t) 9k<ux tk )du 



K tk log x tk , 



(12) 



(13) 
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1 is analogous to a normalized matched filter output for target k 

2 on spatial cell i at time t, and P^iX^) is given in (11) . 

3 Similarly, (7) becomes 

T llz II T U 

t=1 P(x t ) .*t.*\*t-i.k t=1 i=1 

4 (14) 

5 There is an additional term in (14) , but it is omitted here 

6 because it depends on xt,k and not on x tfk , and thus does not 

7 influence the M-step of the EM method. It should be noted at 

8 this point that it is not necessary to have an analytic 

9 expression for S k {v) to utilize (12) and (14) . It is sufficient 

10 to know the values of the set of integrals jj^.^^t^) ' 

11 j = l f --V, for each target Jc. This vector of spectral cell 38 

12 probabilities is a non-parametric description of the target 

13 spectral density sufficient for the problem at hand. 

14 At this stage, specific parametric forms are adopted for 

15 the target and measurement processes. For target k, k=l, . . . , 

16 M, the process evolution is defined by 

17 p_ , l*tt|xt-i f Jt> = N i x tk-> ^t-i^t-i^A-i^) (15) 

a t,Jtpt-l,Jfc 

18 where Af(x; )i,E) is the multivariate normal distribution in x with 

19 mean jx and covariance S. The measurements are characterized by 

12 



2 The covariance matrix R tk relates to the spatial extent, or 

3 spreading, of the energy about its expected location given by 

4 HtjcXtic- Estimates of {n tk \ are obtained using a Lagrange 

5 multiplier technique. The result is 



8 



U 



*tk = £ l. {t) Nte H tk*tk< *t*)du, d7) 

^ 2=1 2 



where 

C7 v 



M 
ic = 0 



S *tki ^ (t) ^( U '' »t*xfcjk. 



= S 2>ij (is) 



9 The last form follows by taking the sum over k innermost and 

10 using Eq. (11) . 

11 Estimates for the state variables are obtained by setting 

12 the gradient of the auxiliary function 0 kx to zero and solving; 

13 however, as in the earlier developments of H-PMHT, an 

14 alternative approach is taken because it exploits the Kalman 

15 filter as an efficient computational algorithm. The details of 

16 the Kalman filter steps are omitted here, however, the synthetic 

17 spatial measurements used in the filter for target k now have 

18 the form 

u , . 

19 *(* = — 2 Vtki J> . « A^ii; H tk x tkf R tk )du, (19) 

2 0 where 
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U / V 

1 v tk = Z ^tJti L rrt M u; H tJc*i*, «dfcjiu- (20) 

2 The synthetic process and measurement noise covariance matrices 

3 used in conjunction with this synthetic measurement are 

4 respectively given by 

p(x ) 

5 Qtk = -ti 0 ^' 0 < t < T — 1 (21) 

Ft+i|| 

6 and 

7 £ tJfc = f tJc , 1 < t < T (22) 

*tk v tk 

8 Let {^Jc} ^ e the set °f estimated mixing proportions and \*tk\ 

9 and |r£jc} define the signal states and width parameters at the /- 

10 th EM iteration. For simplicity and robustness, assume that 

11 {^1=1^} and j^/cj^jt} for a11 t=l # . . .,T in the batch of scans 

12 of the sensor space 35. These restrictions, tantamount to 

13 statistical stationarity , are most often reasonable over the 

14 data intervals of interest. Further, since the spectral density 

15 is never itself required, we will denote the needed integrals by 

16 S k j = js k (v)dv: 

Ej(t) 

17 The method as described below is illustrated in the flow 

18 chart in FIG. 4. At the beginning of the method (the 0-th 

19 iteration) , the mixing proportions j^j^} are initialized so that 

14 



1 n$>0 and + + . . . + n$ = 1 . The signal state sequences- 

2 x$ = jx^, . . . , x^, . . . , x£/cj are initialized with nominal values 

3 for k=l, . . . , M, and the signal widths {r^, . . . , are set 

4 nominally above the expected signal widths so that the tracks 

5 are better able to "see" nearby energy when poorly initialized. 

6 The simple case of x$ = {xq?!' . - - # -x^L' • ■ ■ ' x o*)k\ f (stationary 

7 target) , has proven an effective starting point in many cases. 

8 The process covariance matrices Q t = {o^i/ Qt,2' • ■ • » Qt,M } are 

9 initialized with values tailored to the problem at hand so as to 

10 compromise between smooth tracking and the ability to follow 

11 through aberrant behavior. Typically it is assumed that the 

12 process covariance matrices are constant over time 

13 Q t = Q = {(?!, Q 2 , ■ . . / Q M } • In order to get the iterative estimator 

14 started, initial values are also required for the target state 

15 spectral distributions S = {s 1# S 2 , . . - , S M ) . The simple case of 

16 Sjc = , — —\ has proven an effective starting point for 

[V V V) 

17 estimating the spectra of spatially isolated targets. The above 

18 described initialization of target parameters is step 50 in FIG. 

19 4. 

20 For iterations 7=1,2,..., the following quantities are 

21 computed: 
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1 1. Component spatial cell probabilities for t=l, . . . , T, 

2 i=l f . - . f U, and Jc=0 , 1, . . . ,M: 
3 

4 (23) 
5 

6 2. Component spatial/spectral cell probabilities for 

7 t=l, . . . , T and i = l, . . . , U, j = l, ... V, and k=0 , 1, . . 

8 . ,M: 

9 ifijtft) - *B«t> S kj . (24) 

10 3. Total spatial/spectral cell probabilities for t=l, . . 

11 . , T and i=l, . . . , £/, j'=l, . . . , V: 

12 lg(xy- Z^ifijOfcX (25) 

* = o 

13 4. Total scan probabilities for t=l, . . . , T: 

14 $ = £ ipg(x t ). . (26) 

15 5. Expected sensor space measurement z t ±j for t=l, . , . ,T 

16 i=l, . . . , U, and j = l,. . ., V using equation (4), 
17 

18 6., Spatial cell first moments for t=l, . . . ,T i = l, . . 

19 . , U, and k=l, . . . ,M: 

20 ^ = £ rMr/H^^l&V (27) 
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7. Relative mode contributions for t=l,. . . , T and 
k=0 g 1, . . . ,M: 

^ - I £ w ■ 

i = lj = l p if x V 

8. Estimated mixing , proportions for t=l, . . . , T and 
Jc=0, 1, . . . ,M: 



n tk = 



k>=o* tkl Vt 



9. Synthetic measurements for t=l, . . . , T and &=1, . . 

:</) _ v v z ilf]sAi 

10. Synthetic measurement covariance matrices for t=l, 
• , T and k=l , . . . ,M: 

_ a-" 

11. Synthetic process covariance matrices for t=0,l,. 
. , T-l and k=l, . . . ,M: 

1/) P</)(Jf t> 

where Q tk is treated as a control parameter for the 
process description, and most commonly Q t k=Qk for all 
t=l, . . . ,T in the batch. 



1 12. Estimated spatial states 55 in FIG. 4 

2 = \?oh i *tkr''- < x m\ for t=0,l f . . .,T and k=l t . . 

3 . ,M, using (for computational efficiency) a recursive 

4 Kalman smoothing filter, on the synthetic data z^l with 

5 process and measurement matrices corresponding to 

6 F tkl $1 H tk , 4$ ■ 

7 13. Spatial cell second moments for t=l, . . . , T, i=l, . . 

8 . , U. and k=l , . . . ,M: 



a tki 



£ = ^ (r - H^" 1 ^ N^H^V, 4k X) )^- (33) 



10 14. Average signal width estimates 60 in FIG. 4 for k=l, 

11 • ,M: 
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Z I I ^ tJcJ - (34) 

t = l±=lj = l P ij( X 0 



13 At the completion of iterations / the estimated signal states x (/) 

14 and their width estimates R (/) = {Ri (/) , R 2 (/) , . , . ,Rm (/) } constitute 

15 the track estimate output. 

16 15. Using the track estimate output, compute the average 

17 synthetic spectral power 65 in FIG. 4 for j=l, . . . , V, 

18 and k=l, . . . ,M: 
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2 (Note from (35) that S kj = 1.) 

3 The resulting combination of processed spatial, signal 

4 width and spectral estimates are linked chronologically 70 and 

5 displayed as an image on a computer display screen 75 using 

6 display methods known in the art. 4 

7 The advantages of the present invention over the prior art 

8 are that the resulting method has improved crossing track 

9 performance on sources that have some degree of spectral 

10 distinction. The present invention also avoids the need for 

11 thresholding and peak-picking to produce point measurements. 

12 The spectral estimates (35) may be used to initiate this 

13 estimator when run on subsequent batches of data. In the 

14 preferred embodiment as a new scan is received the oldest scan 

15 of the batch is dropped and the estimation method including the 

16 steps 1 through 14 (formulae (23) through (35)) as stated above 

17 is run in "sliding batch" fashion using the batch length that 

18 provides sufficient smoothing without being unnecessarily long. 

19 . In this case, t in the equations represents the time index 

2 0 within the batch under consideration. The track and spectral 

21 estimates from the previous batch are used as initial values to 

22 start the iterations as outlined. 
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1 The target specific spectral estimates (35) constitute 

2 outputs unto themselves and can be easily computed for arbitrary 

3 track sequences x' and R' used in place of x {/) and R (/) . The 

4 resulting spectral estimates have been termed "track conditioned 

5 spectral estimates," and they serve to give a spectral 

6 characterization to tracks generated via other means. 

7 Obviously many modifications and variations of the present 

8 invention may become apparent in light of the above teachings. 

9 For example :^k(u/x t k) may take a parametric form other than the 

10 normal density given in (16) , g 0 (u) may be other than the uniform 

11 density as implied by (23). While it was shown here that the 

12 spectrum could be handled in a non-parametric form, the methods 

13 are readily extended to treat a parametric spectral description. 

14 In light of the above, it is therefore understood that 

15 within the scope of the appended claims, the invention may be 

16 practiced otherwise than as specifically described. 
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